#ifndef _CENTRALIZED_POLYGON_SYSTEM_H_
#define _CENTRALIZED_POLYGON_SYSTEM_H_

#include <valarray>

#include <boost/shared_ptr.hpp>

#ifdef _CILK_
#include <cilk/cilk.h>
#include <cilk/reducer_opadd.h>
#include <boost/thread.hpp>
#endif

#include <Geex/basics/types.h>
#include "../common/polytope.h"
#include "../common/potential_energy.h"
#include "distance_field_grid_2.h"
#include "discrete_field_system.h" // for polygon sampler
#include "ANN_distance_field.h"

//#define _BOUNDARY_ENERGY_

namespace Field_packing
{
	using namespace Packing_common;
	
	class Centralized_polygon_system : public System_potential_energy
	{
		typedef Packing_polygon::Point_2 Point_2;
		typedef Packing_polygon::Vector	Vector_2;
		typedef CGAL::Line_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Line_2;
		typedef CGAL::Aff_transformation_2<CGAL::Exact_predicates_inexact_constructions_kernel> Transformation_2;	

		const static int GRID_DIM;
		const static int EXPONENT; // 2-4, 4-8, 6-12

	public:
		Centralized_polygon_system(const std::vector<Packing_polygon*>& polygon_set, const Container_box& bndbox, double rm_scale = 1.0, bool up_gravity = false);

		~Centralized_polygon_system() {}

		int optimize(double *trans, bool apply_constraints = false) ;

		double system_energy_at(const double *const x) 
		{
			move_to_new_point(x);
			double me = molecular_energy_at(x);
			//double ge = 0.0;
			//if (g != 0.0)
			double ge = 0.0;
			if (!upward_gravity)
				ge = gravity_energy_at(x); 
			else
				ge = upward_gravity_energy_at(x);

			double be = 0.0;

#ifdef _BOUNDARY_ENERGY_
			if (g != 0.0)
				be = boundary_energy_at(x);
#endif
			return me + ge + be;
			//return me;
		}

		void energy_gradient_at(const double *const x, double *gradient) ;

		// objective function
		double gravity_energy_at(const double *const x);// { return 0.0; }
		double molecular_energy_at(const double *const x) ;

		double upward_gravity_energy_at(const double *const x);

		// debug
		void draw_sample_pnts() const;

	private:
		void move_to_new_point(const double *const x);

		void molecular_energy_gradient_at(const double *const x, double *gradient);
		void gravity_energy_gradient_at(const double *const x, double *gradient);

		void upward_gravity_energy_gradient_at(const double *const x, double *gradient);

		void constraint_at(const double *const x, double *constraint_values) ;
		void constraint_grad_at(const double *const x, double *constraint_grad) ;

		double boundary_energy_at(const double *const x);
		void boundary_energy_gradient_at(const double *const x, double *gradient);

	private:
		
		//Point_2 container_cent;
		double rm; // at which the potential reaches minimum
		bool upward_gravity;
		size_t nb_polygons;

		typedef std::vector<Point_2> Point_set;
		std::vector<Point_set> sample_pnts;
		Point_set centroids;
		std::valarray<double> areas;
		double con_btm; // y coordinate of the bottom of the container
		double con_top;
		std::vector< matrix<double> > prev_trans; // [ cos, sin; tx, ty ]
		std::vector<Point_set> polygon_vtx;
		Point_set bnd_pnts;

		// discretized grid
		Distance_field_grid_2 dfg;

		//kd-tree based field
		//ANN_distance_field_2 adf;

		boost::shared_ptr<Knitro_context> ktr_context;

	};
}
#endif